#if !defined( MATRIX_IO_H)
#define  MATRIX_IO_H
#include <iostream>
#include <vector>
#include <algorithm>
#include <cassert>
#include <stdexcept>
// boost headers
#include <boost/lexical_cast.hpp>
#include <boost/algorithm/string/split_lg.hpp>
#include "is_equal_to.h"
#include "matrix_formats.h"
#include "index_string.h"
#include "matrix.h"
#include "ultoa.h"
#include "progress_indicator.h"
#include "matrix.h"
#include "stlhelper.h"
#include "perlutil.h"
#include "progress_indicator.h"
#include "binary_io.h"
#include "tab_delimited_str.h"

class t_progress_indicator;
namespace bcpp_matrix
{

inline std::string phylip_pad_name(const std::string& name, bcpp_matrix::eMATRIX_IOFORMAT format)
{
	if ((format == bcpp_matrix::eMIO_PHYLIP			|| 
		 format == bcpp_matrix::eMIO_PHYLIP_NULL) &&
		name.length() < 10)
		return (name + std::string(' ', 10 - name.length()));
	else
		return name;
}

template<typename T>
inline std::string get_null_str(T binary_null)
{
	std::ostringstream osstrm;
	osstrm << binary_null;
	return osstrm.str();
}



//________________________________________________________________________________________


// for reading in matrices from stream


//________________________________________________________________________________________



//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	Read and write sparse matrix

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//________________________________________________________________________________________

// writing matrices

//________________________________________________________________________________________
template <typename T> 
void write_matrix(	std::ostream&				os, 
					const t_matrix<T>&			io, 
					eMATRIX_IOFORMAT			format,
					t_progress_indicator& 		dots,
					const T&					null_value = t_matrix_traits<T>::binary_null());
// no dots
template <typename T> 
void write_matrix(	std::ostream&				os, 
					const t_matrix<T>&			io, 
					eMATRIX_IOFORMAT			format,
					const T&					null_value = t_matrix_traits<T>::binary_null())
{
	std::ostringstream ignore;
	t_progress_indicator dots(ignore,  10000000);
	write_matrix(os, io, format,  dots, null_value);
}

//________________________________________________________________________________________

// reading matrices

//________________________________________________________________________________________
// keeping a count of lines
template <typename T>                               		
bool read_matrix(	std::istream&				is, 
					t_matrix<T>&				io, 
					eMATRIX_IOFORMAT			format,
					unsigned&					cnt_lines,
					t_progress_indicator& 		dots,
					const T&					null_value = t_matrix_traits<T>::binary_null());
// no cnt_lines
template <typename T>                               		
bool read_matrix(	std::istream&				is, 
					t_matrix<T>&				io, 
					eMATRIX_IOFORMAT			format,
					t_progress_indicator& 		dots,
					const T&					null_value = t_matrix_traits<T>::binary_null())
{
	unsigned cnt_lines = 0;
	return read_matrix(is, io,  format, cnt_lines, dots, null_value);
}

// no cnt_lines, no dots
template <typename T>                               		
bool read_matrix(	std::istream&				is, 
					t_matrix<T>&				io, 
					eMATRIX_IOFORMAT			format,
					const T&					null_value = t_matrix_traits<T>::binary_null())
{
	std::ostringstream ignore;
	t_progress_indicator dots(ignore,  10000000);
	return read_matrix(is, io, format, dots, null_value);
}





//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	Wrappers to read and write matrix


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

// input stream
template <typename T> 
struct t_matrix_istrm_wrapper
{
	t_matrix<T>&					matrix;
	eMATRIX_IOFORMAT				format;
	unsigned&						cnt_lines;
	t_progress_indicator& 			dots;
	T								null_value;

	t_matrix_istrm_wrapper(	t_matrix<T>&		matrix_,
							eMATRIX_IOFORMAT		format_,
							unsigned&				cnt_lines_,
							t_progress_indicator& 	dots_,
							const T&				null_value_ = 
													t_matrix_traits<T>::binary_null())
		:matrix(matrix_), format(format_), 
			cnt_lines(cnt_lines_), dots(dots_),
			null_value(null_value_){}
};



// output stream
template <typename T> 
struct t_matrix_ostrm_wrapper
{
	const t_matrix<T>&		matrix; 
	eMATRIX_IOFORMAT		format;
	t_progress_indicator& 	dots;
	T						null_value;

	t_matrix_ostrm_wrapper(	const t_matrix<T>&		matrix_,
							eMATRIX_IOFORMAT		format_,
							t_progress_indicator& 	dots_,
							const T&				null_value_ = 
															t_matrix_traits<T>::binary_null())
		:matrix(matrix_), format(format_), dots(dots_),
			null_value(null_value_){}

};


// typedef shortcuts
typedef t_matrix_istrm_wrapper<float> t_matrix_istrm_wrapper_f;
typedef t_matrix_ostrm_wrapper<float> t_matrix_ostrm_wrapper_f;




template <typename T> 
inline std::ostream& operator <<(std::ostream& strm, const t_matrix_ostrm_wrapper<T>& wrapper)
{
	write_matrix(strm, wrapper.matrix, wrapper.format, wrapper.dots, wrapper.null_value);
	return strm;
}
template <typename T> 
inline std::istream& operator >> (std::istream& istrm, t_matrix_istrm_wrapper<T> wrapper)
{
	if (!read_matrix(istrm,  wrapper.matrix,  wrapper.format,
					 wrapper.cnt_lines, wrapper.dots, wrapper.null_value))
		istrm.setstate(std::ios::failbit);
	else
		// in case set at eof
		istrm.clear(istrm.rdstate() & ~(std::istream::failbit));
	return istrm;
}



namespace matrix_io_detail
{
//________________________________________________________________________________________

//	t_id_pair_value

	//	for reading pairs

//________________________________________________________________________________________
template <typename T>
struct t_id_pair_value
{
	T			value;
	unsigned	id_1;
	unsigned	id_2;
	t_id_pair_value(const T& d, unsigned g1, unsigned g2): value(d), id_1(g1), id_2(g2){}
	bool operator < (const t_id_pair_value& o) const
	{
		if (id_1 == o.id_1)
			return id_2 < o.id_2;
		else
			return id_1 < o.id_1;
	}
	t_id_pair_value reversed() const { return t_id_pair_value<T>(value, id_2, id_1);	}
};


//________________________________________________________________________________________

//	process_line_with_pair

//		returns whether diagonal used (id1 == id2)

//________________________________________________________________________________________

template<typename T> 
unsigned process_line_with_pair(const std::string& 					line,
								std::deque<t_id_pair_value<T> >&	id_pair_values,
								t_index_string&						index_string)
{
	std::vector<std::string> data;

	boost::algorithm::split_n( data, line, is_equal_to('\t'), 3);
	//perl_split (line, "\t", perlutil::eIGNORE_EMPTY, data, 3);
	if (data.size() < 3)
		throw std::runtime_error("Reading matrix in pairs format. Three arguments expected per line.\n");
	//
	// names in pair
	//
	unsigned tax1_index = index_string.index(data[0]);
	unsigned tax2_index = index_string.index(data[1]);
	id_pair_values.push_back(t_id_pair_value<T>(t_matrix_traits<T>::from_string(data[2]), tax1_index, tax2_index));
	return (tax1_index == tax2_index ? 1 : 0);
}

//________________________________________________________________________________________

//	pairs_to_matrix

//________________________________________________________________________________________
template<typename T>
void pairs_to_matrix(   const std::deque<t_id_pair_value<T> >&	id_pair_values,
						t_dense_matrix<T>&						matrix,
//						t_progress_indicator&					dots,
						unsigned								matrix_size,
						bool									diagonals_used
						)
{
	bool all_diagonals_missing = true;
	for (typename std::deque<t_id_pair_value<T> >::const_iterator i = id_pair_values.begin(),
		 end = id_pair_values.end(); i != end; ++i)
	{
		unsigned id1 = i->id_1;
		unsigned id2 = i->id_2;
		matrix[id1][id2] =
		matrix[id2][id1] = i->value;
		if (id1 == id2)
			all_diagonals_missing = false;
//		++dots;
	}


	// if all diagonal is missing, assume distance matrix and should be 0
	if (!diagonals_used)
	{
		T zero_value = t_matrix_traits<T>::zero();
		for (unsigned i = 0; i < matrix.size(); ++i)
			matrix[i][i] = zero_value;
	}

}

} //matrix_io_detail

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	matrix in pairs

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
template<typename T>
void read_matrix_from_pairs(std::vector<std::string>::const_iterator	begin,
							std::vector<std::string>::const_iterator	end,
							t_matrix<T>&								data,
							t_progress_indicator&						dots,
							const T&									null_value_flag)
{
	//
	// read pairs into array, with names in hash
	//
	std::deque<matrix_io_detail::t_id_pair_value<T> >	id_pair_values;
	t_index_string										index_names;

	// first string might be matrix name
	if (begin != end && (*begin)[0] == '>')
	{
		data.matrix_name = begin->substr(1);
		begin++;
	}
	unsigned diagonals_used = 0;
	while (begin != end)
	{
		// ignore comments
		if ((*begin)[0] == '#')
		{
			++begin;
			continue;
		}

		++dots;
		diagonals_used += process_line_with_pair(*begin++, id_pair_values, index_names);
	}

	//
	//	make pairs into one big matrix
	//
	data.allocate(index_names.size(), null_value_flag);
	std::copy(index_names.get_deque().begin(), index_names.get_deque().end(), data.names.begin());
	matrix_io_detail::pairs_to_matrix(id_pair_values, data.matrix, 
									  index_names.size(), diagonals_used);

	data.count_missing();
}


template <typename T>
bool read_matrix_from_pairs(	std::istream&				infile,
								t_matrix<T>&				matrix,
								unsigned&					cnt_lines,
								t_progress_indicator&		dots,
								const T&					null_value_flag)
{
	std::deque<matrix_io_detail::t_id_pair_value<T> >	id_pair_values;
	t_index_string				index_names;
	std::string line;
	//dots.set_every(20000);
	bool matrix_read = false;
	unsigned diagonals_used = 0;
	while (infile)
	{
		++cnt_lines;
		++dots;

		//
		//	get matrix name if first line
		//
		if (infile.peek() == '>')
		{
			if (!index_names.size())
			{
				infile.get();
				getline(infile, matrix.matrix_name);
				matrix_read = true;
			}

			//
			//	else matrix name belongs to next matrix
			//
			else
				break;
		}

		// ignore comments
		if (infile.peek() == '#')
		{
			infile.ignore(65536, '\n');
			continue;
		}

		if(!getline(infile, line))
			break;

		if (!line.length())
		{
			// ignore consecutive blank lines
			if (!index_names.size())
				continue;

			// blank lines separate clusters
			else
				break;
		}


		diagonals_used += process_line_with_pair(line, id_pair_values, index_names);
	}

	//std::cerr << id_pair_values.size() << " pairs\n";
	//std::cerr << index_names.size() << " columns\n";
	
	if (index_names.size())
		matrix_read = true;

	if (!matrix_read)
		return false;

	//
	//	make pairs into one big matrix
	//
	matrix.allocate(index_names.size(), null_value_flag);

	std::copy(index_names.get_deque().begin(), index_names.get_deque().end(), 
			  matrix.names.begin());
	matrix_io_detail::pairs_to_matrix(id_pair_values, matrix.matrix, 
									  index_names.size(), diagonals_used);
	matrix.count_missing();

	return true;
}  /* read_matrix_from_stream */



template<typename T>
void write_matrix_to_pairs(std::ostream& os, const t_matrix<T>& dense_matrix, t_progress_indicator& dots)
{
	assert(dense_matrix.size() == dense_matrix.matrix.size());
	assert(dense_matrix.size() == dense_matrix.names.size());

	if (dense_matrix.matrix_name.length())
		os << ">" << dense_matrix.matrix_name << "\n";


	for (unsigned i = 1, sz = dense_matrix.size(); i < sz; ++i)
	{
		for (unsigned j = 0; j < i; ++j)
		{
			if (dense_matrix.matrix[i][j] == dense_matrix.null_value)
				continue;
			++dots;
			os	<< dense_matrix.names[j] << "\t"  
				<< dense_matrix.names[i] << "\t" 
				<< dense_matrix.matrix[i][j] << "\n";
		}
	}
}





//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	matrix in binary

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
template <typename T>
void write_matrix_to_binary(std::ostream& os, const t_matrix<T>& matrix, 
							t_progress_indicator& dots,
							T binary_null)
{
	assert(matrix.size() == matrix.matrix.size());
	assert(matrix.size() == matrix.names.size());


	binary_io::binary_write<unsigned>(os, 12345678);

	binary_io::binary_write(os, matrix.matrix_name);
	binary_io::binary_write(os, matrix.size());


	// use pairwise format if sparse matrix is very sparse
	bool pair_wise_format = matrix.get_count_missing() * 1.0 / 
									matrix.get_count_pairs() > 0.5;
	binary_io::binary_write<unsigned>(os, pair_wise_format ? 1 :0);

	for (unsigned i = 0, sz = matrix.size(); i < sz; ++i)
		binary_io::binary_write(os, matrix.names[i]);


	if (pair_wise_format)
	{
		unsigned cnt_written = 0;
		binary_io::binary_write(os, matrix.get_count_filled());
		// don't assume symmetrical matrix
		for (unsigned i = 0, sz = matrix.size(); i < sz; ++i)
		{
			++dots;
			for (unsigned j = 0; j < sz; ++j)
			{
				++dots;
				if (matrix.is_missing(i, j))
					continue;
				binary_io::binary_write<unsigned short>(os, i);
				binary_io::binary_write<unsigned short>(os, j);
				binary_io::binary_write<T>(os, matrix.matrix[i][j]);
				cnt_written++;
			}
		}
		assert(cnt_written == matrix.get_count_filled());
	}
	else
	{
		// don't assume symmetrical matrix
		for (unsigned i = 0, sz = matrix.size(); i < sz; ++i)
		{
			++dots;
			for (unsigned j = 0; j < sz; ++j)
			{
				++dots;
				if (matrix.is_missing(i, j))
					binary_io::binary_write<T>(os, binary_null);
				else
					binary_io::binary_write<T>(os, matrix.matrix[i][j]);
			}
		}
	}
	// DEBUGG
    //std::cerr << matrix.null_value << "!\n";
}



template <typename T>
bool read_matrix_from_binary(std::istream& 			is, 		
							t_matrix<T>&			matrix,
							unsigned&				cnt_lines,
							t_progress_indicator&	dots,
							const T&				binary_null)
{

	// ignore blank lines at start
	while (is.peek() == '\n')
	{
		++cnt_lines;
		is.ignore(1);
	}
	if (is.eof())
		return false;

	unsigned check = binary_io::binary_read<unsigned>(is);
	if (check != 12345678)
		throw std::runtime_error("Matrix in binary format is corrupted and cannot be read.\n");

	matrix.matrix_name = binary_io::binary_read<std::string>(is);

	unsigned sz = binary_io::binary_read<unsigned>(is);
	matrix.allocate(sz, binary_null);
	bool pair_wise_format = binary_io::binary_read<unsigned>(is) != 0;

	for (unsigned i = 0; i < sz; ++i)
		matrix.names[i] = binary_io::binary_read<std::string>(is);


	if (pair_wise_format)
	{
		unsigned cnt_read = binary_io::binary_read<unsigned>(is);
		while (cnt_read)
		{
			unsigned i = binary_io::binary_read<unsigned short>(is);
			unsigned j = binary_io::binary_read<unsigned short>(is);
			matrix.matrix[i][j] = binary_io::binary_read<T>(is);
			--cnt_read;
			++dots;
		}
	}
	else
	{
		for (unsigned i = 0; i < sz; ++i)
		{
			++dots;
			unsigned j_sz = sz;
			for (unsigned j = 0; j < j_sz; ++j)
			{
				T val = binary_io::binary_read<T>(is);
				if (val != binary_null)
					matrix.matrix[i][j] 		= val;
			}
		}
	}

	matrix.count_missing();
	return true;
}










//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	matrix in square

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
template<typename T>
void write_matrix_to_square(std::ostream&					os, 
							const t_matrix<T>&				matrix, 
							bcpp_matrix::eMATRIX_IOFORMAT	format, 
							t_progress_indicator&			dots,
							T								binary_null)
{
	assert(matrix.size() == matrix.matrix.size());
	assert(matrix.size() == matrix.names.size());


	// no preamble for phylip
	if (format != bcpp_matrix::eMIO_PHYLIP && format != bcpp_matrix::eMIO_PHYLIP_NULL)
	{
		if (matrix.matrix_name.length())
			os << ">" << matrix.matrix_name << "\n";
	}

	os << matrix.size() << "\n";


	std::string str_null = '\t' + get_null_str(binary_null);
	std::string non_null;

	// add null weighting
	if (format == bcpp_matrix::eMIO_PHYLIP_NULL)
	{
		str_null +="\t0";
		non_null +="\t1";
	}

	for (unsigned i = 0, sz = matrix.size(); i < sz; ++i)
	{
		os << phylip_pad_name(matrix.names[i], format);
		// lower doesn't include diagonals
		unsigned j_sz = format == bcpp_matrix::eMIO_LOWER  ? i : sz;
		for (unsigned j = 0; j < j_sz; ++j)
		{
			++dots;
			if (matrix.matrix[i][j] == binary_null)
				os  << str_null;
			else
				os << "\t" << matrix.matrix[i][j] << non_null;
		}
		os << "\n";
	}
}

inline bool do_ignore_comments(std::istream& is, unsigned& cnt_lines)
{
	// ignore comments
	while (is.peek() == '#')
	{
		++cnt_lines;
		is.ignore(65536, '\n');
	}
	if (is.eof())
		return false;

	return true;
}


template <typename T>
bool read_matrix_from_square(	std::istream&					is, 
								t_matrix<T>&					matrix,
								bcpp_matrix::eMATRIX_IOFORMAT	format, 
								unsigned&						cnt_lines,
								t_progress_indicator&			dots,
								const T&						binary_null)
{

	// ignore blank lines at start
	while (is.peek() == '\n')
	{
		++cnt_lines;
		is.ignore(1);
	}
	if (is.eof())
		return false;

	// ignore comments
	if (do_ignore_comments(is, cnt_lines))
		return false;

	bool matrix_read = false;
	if (is.peek() == '>')
	{
		is.ignore(1);
		++cnt_lines;
		getline(is, matrix.matrix_name);
		matrix_read = true;
	}
	if (is.eof())
		return matrix_read;


	// ignore comments
	if (do_ignore_comments(is, cnt_lines))
		return false;


	unsigned sz;
	is >> sz;
	if (!is)
		throw std::runtime_error("Invalid square matrix format.\nThe matrix size on the first line is missing.");
	++cnt_lines;
	matrix.allocate(sz, binary_null);
	std::string str_null = get_null_str(binary_null);

	for (unsigned i = 0; i < sz; ++i)
	{
		// ignore comments
		if (do_ignore_comments(is, cnt_lines))
			return false;

		++dots;
		char endc;
		is >> tab_delimited_str('\n', matrix.names[i], endc);
		++cnt_lines;
		unsigned j_sz = format == bcpp_matrix::eMIO_LOWER  ? i : sz;
		for (unsigned j = 0; j < j_sz; ++j)
		{
			// read first identifier
			std::string value;
			char endc;
			is >> tab_delimited_str('\t', value, endc);
			if (value == "")
				throw std::runtime_error("Invalid square matrix_format.\nEmpty value on line " +
   											ultoa(cnt_lines));

			// phylip replicates
			if (format == bcpp_matrix::eMIO_PHYLIP_NULL)
			{
				std::string replicate_value;
				is >> tab_delimited_str('\t', replicate_value, endc);
				if (replicate_value == "")
					throw std::runtime_error("Invalid square matrix_format.\nEmpty replicate number on line " +
												ultoa(cnt_lines));
				int cnt_replicates = t_matrix_traits<int>::from_string(replicate_value);
				if (!cnt_replicates)
					continue;
			}


			if (value != str_null)
				matrix.matrix[i][j] = t_matrix_traits<T>::from_string(value);

			// lower format is the only one which is symmetrical
			if (format == bcpp_matrix::eMIO_LOWER)
				matrix.matrix[j][i] =matrix.matrix[i][j];
			
		}
		if (!is)
			throw std::runtime_error("Invalid square matrix format.");
	}

	// lower format has zeroed diagonals
	T zero_value = static_cast<T>(0U);
	if (format == bcpp_matrix::eMIO_LOWER)
		for (unsigned i = 0; i < matrix.matrix.size(); ++i)
			matrix.matrix[i][i] = zero_value;

	matrix.count_missing();
	return true;
}














//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	all formats

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
template<typename T>
void write_matrix(	std::ostream&					os, 
					const t_matrix<T>&				matrix,
					eMATRIX_IOFORMAT				format,
					t_progress_indicator& 			dots,
					const T&						binary_null)
{
	switch (format)
	{
		case eMIO_BINARY:
			write_matrix_to_binary(os, matrix, dots, binary_null);
			break;
		case eMIO_PAIRS :
			write_matrix_to_pairs(os, matrix, dots);
			break;
		case eMIO_LOWER :
		case eMIO_SQUARE:
		case eMIO_PHYLIP:
		case eMIO_PHYLIP_NULL:
			write_matrix_to_square(os, matrix, format, dots, binary_null);
			break;
		default:
			assert(1);
	}
}

template <typename T>
bool read_matrix(	std::istream&				is,
					t_matrix<T>&				matrix,
					eMATRIX_IOFORMAT			format,
					unsigned&					cnt_lines,
					t_progress_indicator& 		dots,
					const T&					binary_null)
{
	try
	{
		switch (format)
		{
			case eMIO_BINARY:
				return read_matrix_from_binary(is, matrix, cnt_lines, dots, binary_null);
			case eMIO_PAIRS :
				return read_matrix_from_pairs(is, matrix, cnt_lines, dots, binary_null);
			case eMIO_LOWER :
				return read_matrix_from_square(is, matrix, format, cnt_lines, dots, binary_null);
			case eMIO_SQUARE:
			case eMIO_PHYLIP:
			case eMIO_PHYLIP_NULL:
				return read_matrix_from_square(is, matrix, format, cnt_lines, dots, binary_null);
			default:
				assert(1);
		}
	}
	catch (std::runtime_error& error)
	{
		throw std::runtime_error(std::string(error.what()) + "\nThe error occurred at line # " + ultoa(cnt_lines + 1));
	}

	throw std::logic_error("Write matrix: unknown format\n");
	return true;
}


};




#endif	//MATRIX_IO_H

